%--------------------------------------------------------------------------
% computes the time-dependent eigenvalues numerically and using the method
% of multipe scales. that is, it computes alpha(t) where alpha is from
%
%              c = c(t=0) * exp(int_0^t alpha(u) du)
%
%--------------------------------------------------------------------------

function [t, alpha_as, A] = time_eigs(k, p)
close all;

if nargin == 1
    p = params;
end

lin_soln = stab(k, p);

[c, dcdt] = deval(lin_soln, lin_soln.x, p.N);
alpha = dcdt ./ c;

t = lin_soln.x';


alpha_as = small_k(t, k, p);
% alpha_as = Ma_delta0(t, k, p);
A = c(1) * exp(cumtrapz(t, alpha_as));

if nargin == 1
    
    subplot(2,1,1);
    plot(t, c, 'k', t, A, 'r*', 'linewidth',1);
    xlabel('$t$', 'interpreter','latex','fontsize',12);
    ylabel('$\hat{c}(h(t),t)$', 'interpreter','latex','fontsize',12);
    l = legend('numerical','asymptotic','orientation','horizontal','location','best');
    set(l, 'interpreter','latex','fontsize', 11);
    
    subplot(2,1,2);
    
    if (all(alpha > 0) || all(alpha < 0))
        %     semilogy(t, alpha, 'k',t, ev*ones(size(t)), 'r--',  t, alpha_as, 'r*');
        semilogy(t, alpha, 'k', t, alpha_as, 'r*');
        
    else
        %     plot(t, ev*ones(size(t)), 'k--', t, alpha, 'k', t, alpha_as, 'r*');
        plot(t, alpha, 'k', t, alpha_as, 'r*', 'linewidth',1);
        
    end
    xlabel('$t$', 'interpreter','latex','fontsize',12);
    ylabel('$\alpha$','interpreter','latex','fontsize',12);
    l = legend('numerical', 'asymptotic','orientation','horizontal','location','best');
    set(l, 'interpreter','latex','fontsize', 11);
    
    % ylim([-0.1 0.1]);
end

%--------------------------------------------------------------------------
% alpha for Ma = O(1)

function alpha_as = Ma_delta0(t, k, p)

f0 = -(1 + lambertw(-p.beta / (p.beta - 1) * exp(((-p.beta + p.delta * t) / (p.beta - 1))))) * (p.beta - 1);
A0 = (1 - p.beta) * (1 - 1 ./ f0);

if (k < sqrt(p.delta))
    w = @(z,i) p.Ma * z.^2 .* (k^2 / 4 * (1 - z / f0(i)) + p.delta / p.Ma / 2 / f0(i)^2 * (3 - z / f0(i)));
else
    A = w_coeffs(1, f0, k, p);
    w = @(z,i) (A(i,1) * z + A(i,2)) .* cosh(k*z) + (A(i,3) * z + A(i,4)) .* sinh(k*z);
end

dv1_dz = @(z,i) -(1 - p.beta - A0(i)) * (p.beta + A0(i)) * z / f0(i);

for i = 1:length(t)
    g0(i) = 1 / f0(i) * quad(@(z) w(z, i) .* dv1_dz(z, i), 0, f0(i));
end

alpha_as = -k^2 + p.delta * (-g0' - (1 - 2 * (p.beta + A0)) ./ f0);

%--------------------------------------------------------------------------
% alpha for k = O(1 / sqrt(Ma))

function alpha_as = small_k(t, k, p)

f0 = -(1 + lambertw(-p.beta / (p.beta - 1) * exp(((-p.beta + p.delta * t) / (p.beta - 1))))) * (p.beta - 1);
A0 = (1 - p.beta) * (1 - 1 ./ f0);

alpha0 = -k^2;
alpha1 = (1 - p.beta - A0) .* (p.beta + A0) .* p.Ma .* k^2 .* f0.^2 / 80 - (1 - 2 * (p.beta + A0)) ./ f0;
alpha2 = 11 / 40 * (1 - p.beta - A0) .* (p.beta + A0);

tmp = 1 + (f0 - 1) / p.beta;
alpha1 = p.Lambda / 80 * k^2 * tmp - (1 - 2 * p.beta * tmp ./ f0) ./ f0;
alpha2 = 0;

alpha_as = alpha0 + alpha1 * p.delta + 0*alpha2 * p.delta^2;